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important  near  the  sharp  leading  edges,  are  simulated  by  a 
Kutta  condition  applied  at  the  leading  edges  of  the  wings. 

The  rest  of  the  flow  field  is  essentially  controlled  by  the 
inviscid  equations.  The  equations,  written  in  conservation 
form  in  generalized  curvilinear  coordinates,  are  approximated 
by  MacCormack's  second-order  accurate  predictor-corrector 
algorithm.  The  flow  tangency  conditions  at  the  body  surface 
are  satisfied  by  Abbett's  scheme  and  the  outer  bow-shock 
position  by  the  Rankine-Hugoniot  jump  relations.  Internal 
shock  waves  or  tangentially  discontinuities  are  captured. 

The  slip  surface  emanating  from  the  leading  edge  of  the 
sharp  wing  excited  nonlinear  instabilities  in  the  MacCormack's 
scheme,  which  were  stabilized  by  special  flow  dependent 
fourth-order  damping  terms. 

Comparisons  have  been  made  between  predicted  pressure 
distributions  and  measured  pressure  distribution  for  an  ®  =  1 
delta  wing  at  a  =  10°  at  =  3.0.  Good  agreement  is  obtained. 
A  solution  was  also  obtained  for  the  same  wing  mounted  on  a 
body  of  revolution  at  the  same  angle  of  attack  and  Mach  number. 
Large  losses  of  favorable  wing-body  interference  were  predicted 
which  are  in  good  agreement  with  experiment. 
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SUMMARY 


This  paper  addresses  the  problem  of  computing  the  steady 
inviscid  supersonic  flows  about  thin  wings  having  sharp  sub¬ 
sonic  leading  edges  and  leading-edge  separation.  For  such 
wings  the  flow  on  the  suction  side  tends  to  spiral  into  a  vor¬ 
tex.  Both  wings  alone  or  wings  in  the  presence  of  the  body  are 
considered. 

Previous  analytical  studies  to  solve  these  flow  fields 
have  used  the  leading-edge  suction  analogy,  linear  slender-wing 
theory,  or  detached  flow  methods.  These  studies  are  basically 
inviscid.  Viscous  effects  have  been  included  by  Vigneron, 
et  al.  in  a  numerical  study  of  supersonic  flows  over  delta 
wings  with  sharp  subsonic  leading  edges  using  the  Navier- 
Stokes  equations.  To  obtain  a  more  efficient  procedure  than 
Vigneron' s  et  al.,  we  use  the  steady  Euler  equations  as  the 
basic  governing  equations  as  opposed  to  the  more  complicated 
Navier-Stokes  equation.  The  viscous  effects,  important  near 
the  sharp  leading  edges,  are  simulated  by  a  Kutta  condition 
applied  at  the  leading  edges  of  the  wings.  The  rest  of  the 
flow  field  is  essentially  controlled  by  the  inviscid  equations. 
The  equations  are  written  in  conservation  form  in  generalized 
curvilinear  coordinates.  The  equations  are  approximated  by 
MacCormack's  second-order  accurate  predictor-corrector 
algorithm.  The  flow  tangency  conditions  at  the  body  surface 
are  satisfied  by  Abbett's  scheme  and  the  outer  bow-shock 
position  by  the  Rankine-Hugoniot  jump  relations.  Any  internal 
shock  waves  or  tangentially  discontinuities  are  captured  by  the 
scheme.  It  was  found  that  the  slip  surface  emanating  from  the 
leading  edge  of  the  sharp  wing  excited  nonlinear  instabilities 
in  the  MacCormack's  scheme.  The  addition  of  special  flow 
dependent  fourth-order  damping  terms  stabilized  the  scheme. 
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Comparisons  have  been  made  between  predicted  pressure 
distributions  and  measured  pressure  distribution  for  an  M =  1 
delta  wing  at  o  =  10°  at  MB  =  3.0.  Good  agreement  is  obtained. 
A  solution  was  also  obtained  for  the  same  wing  mounted  on  a 
body  of  revolution  at  the  same  angle  of  attack  and  Mach  number. 
Large  losses  of  favorable  wing-body  interference  were  predicted 
which  are  in  good  agreement  with  experiment. 
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1.  INTRODUCTION  AND  BACKGROUND 


As  a  result  of  past  work  (ref.  1)  done  at  Nielsen 
Engineering  &  Research,  Inc.  (NEAR)  for  the  Office  of  Naval 
Research  (ONR)  to  produce  engineering  prediction  methods  for 
missiles  up  to  angle  of  attack  of  50°  and  Mach  numbers  up  to 
3.0,  several  important  problems  in  high  angle  of  attack  aero¬ 
dynamics  have  emerged.  It  is  the  purpose  of  this  report  to 
describe  these  problems  and  to  attempt  their  solution  using 
the  Euler  equations.  While  these  problems  have  arisen  in 
connection  with  missile  aerodynamics,  they  are  equally  impor¬ 
tant  for  airplanes. 

The  present  work  is  a  logical  extension  of  the  work  NEAR 
has  carried  out  for  ONR  over  the  past  four  years.  During  the 
first  two  years,  paneling  (inviscid)  methods  together  with 
vortex  theory  were  used  to  develop  pressure  predictive  tech¬ 
niques  for  complete  missiles  at  supersonic  speeds  and  angles 
of  attack  to  20°.  In  the  third  and  fourth  years  engineering 
methods  were  developed  to  predict  forces  and  moments  acting 
on  canard  cruciform  missile  to  50°  angle  of  attack.  These 
methods,  based  on  data  base  and  rational  modeling,  produced 
useful  methods  of  reasonable  accuracy,  but  also  uncovered 
important  aerodynamic  effects  which  need  more  study  in  their 
own  right.  In  order  to  extend  the  range  of  the  engineering 
design  codes  and  to  improve  the  rational  modeling  approach, 
it  is  necessary  to  undertake  further  work.  It  is  believed 
that  such  work  should  utilize  recent  advances  in  viscous  and 
inviscid  numerical  computation  techniques  supplemented  by 
further  experimental  work. 

Two  of  the  major  problems  which  have  been  uncovered  are 
the  adverse  effect  of  wing-body  interference  on  wing  lift  at 
high  angles  of  attack  and  the  special  behavior  of  body  vortices 
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at  supercritical  crossflow  Mach  numbers.  This  report  addresses 
the  first  problem  only. 

On  a  wing-body  combination  at  subsonic  and  moderately 
supersonic  speed,  the  effect  of  the  body  on  the  wing  is  to 
increase  the  wing  lift  so  long  as  the  angle  of  attack  is  small. 
As  an  example,  consider  the  upwash  in  a  crossflow  plane  normal 
to  the  body.  As  it  goes  around  the  body,  the  air  speeds  up  and 
produces  an  angle  of  attack  equal  to  2a  at  the  side  of  the  body 


Putting  the  wing  panels  in  this  field  causes  an  increased  lift 
over  what  they  would  develop  if  the  panels  were  made  into  a 
wing  at  angle  of  attack  a.  The  wing-body  interference  factor 
is  a  measure  of  this  interference. 

Lift  of  wing  panel  attached  to  fuselage  (at  a) 

HW  Lift  of  wing  panel  in  wing  alone  (at  a) 

According  to  potential  flow  theory,  runs  from  1  for  no 
fuselage  up  to  2  for  a  small  fin  on  a  large  fuselage  (ref.  2). 

In  the  previous  work,  we  have  determined  Kw  from  experi¬ 
mental  data  for  a  series  of  wings  and  Mach  numbers  up  to  45° 
angle  of  attack.  Figure  1  shows  these  experimental  results 
for  a  delta  fin  of  aspect  ratio  1  mounted  on  a  body  whose 


diameter  is  one-half  the  total  span  of  the  wing-body.  As  the 
angle  of  attack  is  increased,  the  value  of  decreases  sharply 
indicating  less  favorable  interference  of  the  body  on  the  wing. 
At  high  Mach  numbers  the  adverse  interference  is  substantially 
increased. 


It  is  of  interest  to  examine  the  environment  in  which  a 
fin  or  wing  might  operate  at  high  angles  of  attack.  We  are  in 
a  position  to  calculate  the  inviscid  flow  in  which  such  a  fin 
or  wing  might  operate  using  a  computer  program  (refs.  3  and  4) 
for  solving  the  Euler  equations  developed  at  NASA/Ames  Research 
Center.  Figure  2  shows  such  a  flow  for  a  supersonic  crossflow 
Mach  number  with  the  computer  program  operating  in  the  shock¬ 
capturing  mode.  Within  the  bow  shock  from  the  nose  of  the  body 
the  flow  in  the  crossflow  plane  is  mostly  subsonic.  However 
there  are  embedded  supersonic  flow  regions  bounded  by  M  =  1 
lines  and  crossflow  shocks.  The  flow  as  shown  represents  an 
inviscid  flow  with  rotation  caused  by  the  bow  shock  curvature. 
The  vortical  singularity  line  exhibited  by  the  flow  is  charac¬ 
teristic  of  such  solutions.  In  the  real  fluid  case,  the  wake 
of  the  body  will  modify  the  flow  and  may  eliminate  the 
"vortical  singularity". 

It  is  clear  that  inserting  a  fin  in  the  body  flow  which 
is  shown  in  figure  2  can  cause  significant  changes  in  the  fin 
aerodynamic  force  and  moment  over  what  they  would  be  for  low- 
speed  crossflow  past  the  body.  In  fact  this  altered  flow  field 
accounts  in  the  main  for  the  considerable  reduction  in  for 
this  condition  as  shown  by  figure  1. 

This  report  addresses  the  problem  of  computing  the  steady 
inviscid  supersonic  flows  about  wings  having  subsonic  leading 
edges  with  leading-edge  separation.  For  such  wings  the  flow  on 
the  suction  side  tends  to  spiral  into  a  vortex  above  the  wing. 
This  vortex  provides  for  lift  augmentation  at  low  supersonic 
speeds.  Previous  analytical  studies  to  solve  these  flow  fields 


have  used  the  leading-edge  suction  analogy,  linear  slender-wing 
theory,  or  detached  flow  methods  {ref.  5).  These  studies  are 
basically  inviscid.  Viscous  effects  have  been  included  by 
Vigneron,  et  al.  (ref.  6)  in  a  numerical  study  of  supersonic 
flows  over  delta  wings  with  sharp  subsonic  leading  edges  using 
the  Navier-Stokes  equations . 

Here  we  wish  to  calculate  flows  around  wings  and  wing-body 
combinations  with  sharp  leading  edges.  To  obtain  a  more  effi¬ 
cient  procedure  than  Vigneron' s  et  al.,  we  use  the  inviscid 
steady  Euler  equations  as  the  basic  governing  equations  as 
opposed  to  the  more  complicated  Navier-Stokes  equation.  The 
viscous  effects,  important  near  the  sharp  leading  edges,  are 
simulated  by  a  Kutta  condition  applied  at  the  leading  edges  of 
the  wings.  The  rest  of  the  flow  field  is  essentially  controlled 
by  the  inviscid  equations.  To  further  simplify  the  calculations, 
we  consider  cases  for  which  the  axial  Mach  number  is  supersonic. 

A  brief  outline  of  the  general  approach  of  studying  the 
wing-body  interference  problem  for  steady  supersonic  flows  is 
now  given. 

The  steady  Euler  equations  are  written  in  conservation 
form  in  generalized  curvilinear  coordinates.  The  general  coor¬ 
dinate  system  is  fitted  between  the  body  and  the  outer  bow 
shock.  The  coordinate  system  is  fitted  to  the  body  by  a  con¬ 
formal  transformation  which  transforms  the  wing-body  or  wing 
alone  to  a  unit  circle.  To  fit  the  outer  shock  to  the  coordi¬ 
nate  system,  the  radial  distance  is  normalized  with  respect  to 
the  distance  between  the  outer  shock  and  the  body.  A  typical 
mesh  is  shown  in  figures  3(a)  and  (b)  for  the  planar  delta  wing. 
The  equations  are  approximated  by  MacCormack's  second-order 
accurate  predictor-corrector  algorithm  (refs.  3  and  4) .  The 
flow  tangency  conditions  at  the  body  surface  are  satisfied  by 
Abbett's  scheme  (ref.  7)  and  the  outer  bow- shock  position  by 
the  Rankine-Hugoniot  jump  relations.  Any  internal  shock  waves 
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or  tangentially  discontinuities  are  captured  by  the  scheme.  It 
was  found  that  the  slip  surface  emanating  from  the  leading  edge 
of  the  sharp  wing  excited  nonlinear  instabilities  in  the 
MacCormack's  scheme.  The  addition  of  special  flow  dependent 
fourth-order  damping  terms  stabilized  the  scheme. 

Because  of  the  singularity  of  the  conformal  transformation, 
the  transformed  Euler  equations  become  indeterminate  as  the  wing 
tip  is  approached.  A  separate  specification  of  the  primitive 
flow  variables  at  the  wing  tip  is  required.  This  is  done  by 
satisfying  the  Kutta  condition,  that  is,  the  flow  at  the  wing 
tip  is  tangent  to  the  wing  surface.  The  pressure,  entropy,  and 
flow  direction  (in  the  plane  of  the  wing)  are  obtained  by 
averaging  the  values  obtained  numerically  just  above  and  below 
the  wing  tip.  Constant  total  enthalpy  completes  the  specifica¬ 
tion  of  the  wing  tip  flow  variables.  For  the  wing-alone 
results,  no  wing  thickness  is  used.  The  singularity  of  the 
transformation  at  the  wing-body  juncture  is  avoided  by  moving 
a  distance  e  off  the  unit  circle. 

In  the  following  four  sections,  the  governing  transformed 
equations,  the  body  geometry,  the  numerical  scheme  and  the 
boundary  and  initial  conditions  are  described.  The  numerical 
results  are  presented  and  discussed  in  the  sixth  section.  In 
the  seventh  section  conclusions  and  recommendations  are  made. 

The  final  section,  an  appendix,  describes  briefly  the  graphics 
required  to  obtain  the  crossflow  particle  trajectories  and  the 
conical  streamlines. 

2.  STEADY  EULER  EQUATIONS 

The  conservation-law  form  of  the  fluid  dynamic  equation 
for  steady,  inviscid,  three-dimensional  compressible  flow  of 
an  ideal  gas  (steady  Euler  equations)  in  Cartesian  form  are 
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where  E,  P ,  and  6  are  defined  as 


pu 
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pu  +  kp 
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puv 

A 

puw 
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2 
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puw 

pvw 

2 

pw  +  pk 

—  _ 

Equation  (1)  represents  the  conservation  of  mass  and 
momentum.  The  pressure  and  density  are  normalized  with  respect 
to  the  stagnation  conditions  and  the  Cartesian  velocity  com¬ 
ponents  (u,v,w)  with  respect  to  the  maximum  adiabatic  velocity 
where  k  =  2y/(Y-l)  and  Y  is  the  ratio  of  the  specific  heats. 

The  system  of  equations  is  closed  by  the  integrated  form  of 
the  steady  energy  equation  which  in  nondimensional  form  is 

p  =  p(l  -  u2  -  v2  -  w2)  (2) 

For  a  given  free-stream  Mach  number  and  angle  of  attack  a, 
the  remaining  free-stream  variables  in  nondimensional  form  are 
given  by 

?  -Y/(Y-1) 
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-  -1/(Y-1) 

P„  =  {1  +  [(Y-1)/2]M* 

“=o  =  0 

=  q  sin  a 

w  =  q  cos  a 
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q  =  (u2  +  v2  +  w2)1/2 
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To  obtain  a  surface-oriented  coordinate  system,  the  system 
(1)  is  transformed  from  the  Cartesian  space  (x,y,z)  into  another 
arbitrary  curvilinear  system  (z,r,<j>)  where  r  =  1  defines  the 
body  surface.  The  general  transformation  is  (the  particular 
transformations  considered  are  given  in  section  3) 


z  =  z 

r  =  r (x,y, z)  (3) 

4>  =  <f>  (x,y ,  z) 

The  transformed  equation,  obtained  by  Viviand's  (ref.  8) 
procedure,  are  in  slightly  rearranged  form 

+  +  (4) 

where 
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U  =  w 

V  =  ru  +  rv  +  rw 

x  y  z 

W=0u+<tv+<bw 


(5) 


where  U,  V,  and  W  are  the  contravariant  velocities  written 
without  the  metric  normalization.  The  metric  terms  are  obtained 
from  the  chain  rule  expansion  of  xr,  yr,  etc.  and  solved  for 

rx'  ry'  etc*  to  give 
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(6) 


and 
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J  -  V*  -  Vr 


where  J  is  the  Jacobian  of  the  transformation  from  the  arbitrary 
curvilinear  space  (z,r,<|>)  to  the  Cartesian  space  (x,y,z), 
figure  4.  It  should  be  pointed  out  that  both  coordinate 
systems  are  left-handed,  the  reason  being  that  the  parent 
code  (refs.  3  and  4)  used  for  this  study  was  written  for  a 
left-handed  system. 

To  fit  the  outer  bow  shock  wave,  the  outer  mesh  boundary 
must  coincide  with  the  bow  shock.  Since  this  bow  shock  is  a 
variable  three-dimensional  surface  it  is  necessary  to  introduce 
another  transformation  which  normalizes  the  distance  between 
the  body  boundary  and  the  how  shock  surface.  The  location  of 
the  bow  shock  surface  is  determined  by  the  Rankine-Hugoniot 
conditions  during  the  course  of  the  numerical  computation 
described  in  a  later  section.  At  the  same  time  it  is  desirable 
to  have  arbitrary  clustering  functions  in  the  transformation  so 
that  mesh  points  can  be  concentrated  near  the  body  surface, 
wing  tip,  or  wing-body  juncture  for  increased  resolution  in 
areas  of  rapid  changes  of  the  flow  variables  or  the  previously 
mentioned  transformation  metrics. 
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p . . pg*p«Pl^pia«piiiiinqBp 


The  equations  of  the  independent  variable  transformation 


are 


5  =  2 

t  =  h(5> 
n  =  f  <  ) 


(7) 


where  h  and  f  are  clustering  functions 


C  -  (r  -  rb)/(rs  -  rb> 


and 


rb  ■  rb(2,4»),  the  body  surface  radius 


r  =  r _(z,4>),  the  outer  shock  radius 
s  s 


to  be  determined  as  part  of  the  numerical  solution  procedure. 
The  derivatives  are 


A  =  A  +  T  $  A 

3z  35  52  3t 


A  =  x  5  A 

3r  r  9t 


(8) 


_3_  _  _  r  A  +  n  A 
3<}>  “  51  3t  3r) 

The  system  (1)  is  now  in  weakly  conservative  form  (i.e.  no 
longer  strongly  conservative  in  the  spirit  of  Viviand) 


£e  +  1Fp  Ag  +  h=  0 


(9) 


where 


E  =  E 


F  =  +  i-rF  +  ^G} 

g  =  n#Q 
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i| 

i 


and 


We  also  have 


T 


5 


’55 


d2h(£) 


*  1/(rs  * 


V 


KZ  =  “{rbz  +  S<r«7  -  rK7)}  •  5 


sz  bz 


(10) 


Ct  ’  ’{rb,  +  «<r.*  -  rb*))  •  t, 


„  =  and  ,  £im 

<p  d$  '<fxp  2 

d<p 


The  functions  h  and  f  are  the  clustering  transformations 
in  the  r  and  <p  directions,  respectively.  The  normalization 
between  the  body  surface  and  the  shock  is  given  by  the  variable 
5. 

The  equation  is  no  longer  in  strongly  conservative  form 
(i.e.  H  =  0)  to  simplify  the  decoding  of  the  dependent  E  vector. 
The  finite  difference  form  of  eq.  (9)  is  integrated  with  re¬ 
spect  to  the  hyperbolic  coordinate  c  to  yield  values  of  E. 

The  physical  flow  variables  p,  p,  u,  v,  w  must  be  decoded  from 
the  components  e^  of  E.  Explicit  expressions  are  possible  for 
the  physical  flow  variables  if  the  system  is  weakly  conservative. 


t 

I 


if 


!  ■ 

Ij 


t 

i 


II 


£ 


I 


i 
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This  is  not  possible  in  the  strongly  conservative  form  since 
the  E  vector  would  have  a  form  similar  to  F  in  eq.  (9) . 


The  physical  flow  variables  are  obtained  by  the  solution 
of  five  simultaneous  nonlinear  equations  consisting  of  the  four 
components  ei  and  the  integrated  energy  equation  (eq.  (2)). 

The  velocity  components  u  and  v  are  given  by 


■  e4/el 
=  e3/e 


The  ei  along  with  eqs.  (11)  and  the  Jacobian  are  used  to 
obtain  an  implicit  relation  in  w  from  eq.  (2) . 


(l-k)w2 


w  + 


(12) 


The  decoding  procedure  is  reduced  to  finding  the  roots  of 
the  quadratic  eq.  (12) .  Two  roots  exists  corresponding  to  a 
subsonic  and  a  supersonic  flow  in  the  z— direction.  The  root 
corresponding  to  supersonic  flow  is  the  desired  solution. 


w  = 


2e. 


1 

(l-k) 


( 


e~  +  {e,  -  4(l-k)k(ef  - 


e,  - 


e4>> 


1/2 


(13) 


The  density  is  obtained  from 

p  =  e^  •  J/w  (14) 

and  the  pressure  from  eq.  (2). 

This  decoding  procedure  gives  explicit  expressions  for 
the  physical  flow  variables.  For  the  strongly  conservative 
form  of  the  governing  equations  a  Newton-Raphson  procedure  is 
necessary  to  determine  the  flow  variables.  The  additional 
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computing  time  was  not  thought  to  be  worthwhile  simply  to 
maintain  the  strongly  conservative  form.  Apparently  the  cor¬ 
rect  jump  conditions  are  still  obtained  with  the  equations  in 
weakly  conservative  form  (ref.  3). 


3 .  BODY  GEOMETRY 


In  this  study  two  types  of  body  geometries  are  considered: 
a  planar  wing  alone  and  a  planar  mid-wing  on  a  cylindrical  body. 
The  latter  configuration  is  shown  in  figure  4.  The  wing  lies 
in  the  x-z  plane.  Both  the  body  radius  and  wing  semispan  are 
arbitrary  functions  of  the  axial  coordinate  z.  By  a  conformal 
transformation  the  wing-body  configuration  is  transformed  into 
a  cylinder  of  unit  radius.  The  transformation  is  given  more 
conveniently  in  the  terms  of  the  transformed  space  variables 
(z,r,<J>).  The  mapping  between  the  transformed  space  and  the 
physical  space  (x,y,z)  is  given  by 

g  =  |jz  +  i±  /(Z  +  |)2  -  16B2/R2  )  (15) 

z  =  z 


where 


g  =  x  +  iy 


i 

R 


Z 


/=! 


Rw(z)  =  B2(z)/Rw(z) 

i  ( d>— tt/2 )  i6 

re  y  '  =  re 


R  >  B 

W  — 


and  the  positive  (negative)  sign  in  front  of  radical  applies 
to  the  upper  (lower)  half  plane. 

This  single  transformation  covers  all  the  configurations 
considered  in  this  report.  Circular  bodies  alone,  planar  wings 
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alone,  and  planar  mid-wing,  circular-body  combinations  are  ob¬ 
tained  simply  by  setting  Rw=B,  B=0,  and  >  B,  respectively. 

The  transformation  metrics  are  easily  obtainable  by  dif¬ 
ferentiation  of  the  above  mapping. 


Separating  real  and  imaginary  parts  of  gz,  gr,  and  g^  obtains 

xz,yz,  xr,yr,  and  x^,y^,  respectively 

Substituting  these  into  eq.  (6)  obtains  the  requisite  metrics 
for  eq.  (4) . 

The  transformation  is  well  behaved  (i.e.  J  1  /  0  or  ») 
everywhere  except  at  the  wing  tip  and  at  the  wing-body  juncture. 
At  the  wing  tip  all  the  metrics  and  the  Jacobian  vanish,  leav¬ 
ing  eq.  (4)  indeterminate.  This  requires  special  procedures  at 
the  wing  tip  as  discussed  in  section  5.  At  the  wing-body 
juncture,  the  inverse  of  the  radical  of  eqs.  (16)  becomes 
infinite.  This  singular  behavior  is  due  to  the  abrupt  jump  of 
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the  surface  slope  between  the  body  and  wing.  If  the  analytic 
expressions  (eq.  (It))  for  the  metrics  are  replaced  by  their 
finite  difference  equivalent,  the  singularity  is  automatically 
avoided  and  no  special  treatment  is  required  at  these  points. 
However  a  different  procedure  is  used  here  as  described  in 
section  5.E. 

4.  FINITE  DIFFERENCE  EQUATION 

The  flow  equations,  eq.  (9),  are  approximated  by 
MacCormack's  second-order,  predictor-corrector  scheme.  Since 
the  outer  bow  shock  wave  is  fitted  and  the  internal  shock  waves 
or  tangential  discontinuities  are  captured,  it  is  essential 
that  an  efficient  numerical  scheme  be  used.  Kutler  has  found 
(ref.  9)  MacCormack's  scheme  to  be  one  of  the  most  efficient 
and  easily  implementable  schemes  available. 

The  numerical  algorithm  can  be  written  for  the  field 
points  (4  i  j  i  NT2-1;  3  ^  k  i  NPHI,  see  fig.  5)  as 


2n+l  s  En  -  ^  (Fn  -  Fn  )  -  ff.n  n  . 

3  »k  j,k  At  lFj+l,k  Fj,k)  An  (Gj,KPL  Gj,KPR} 


A£n+V?  .  +  AEn  .  I  +  AE1?  .  I 
3,k  D,k|m  J>k|r 


for  the  predictor  step  and  as 


rn+l  _  1  f  n  n+1  _ 

'j.k  ‘  2  (_Ej.k  +  Ej,k 


IFn+l  .  Fn+1  , 

At  j  ,k  j-l,k; 


_n+l  _n+l  .  ■  _n+l. 


.n+1 1 


[IJ 


for  the  corrector  step,  where 


KPL  =  k  +  NFLIP 
KPR  =  k  -  1  +  NFLIP 
KCL  =  k  +  1  -  NFLIP 
KCR  =  k  -  NFLIP 

E”  =  E  Un,  jAi ,  kAn) 

J  t  K 

Fj#k  =  F<Ejfk'  c"»  kArl) 
p"*£  =  F(e”+J,  ;n+Acn+1,  j At,  kAn) 

and 

Cn  =  c  +  l  hr.1 
°  £=1 

The  value  of  NFLIP  alternates  cyclically  between  one  and 

zero  with  the  integration  steps  to  obtain  unbiased  results. 

The  increments  At  and  An  are  the  mesh  spacings  in  the  radial 

and  meridional  directions  and,  A£n  is  the  marching  step  size 

between  the  n-1  and  nth  computational  planes.  The  terms 

AE1?  .  I  and  AE*?+f  |  are  fourth-order  damping  terms  intro- 
j  ,  k  |m,  r  ] ,k|m,r 

duced  to  control  the  severe  nonlinear  instabilities  experienced 
near  the  leading  edge  of  the  wing .  The  damping  terms  are  given 


i 


where 


and 


c 

m 

pj,k+£+l 

-  2p3 

,k +1  + 

P j , k+£-l 

2 

p j ,k+£+l 

+  2  p 

J 

,k +Z  + 

pj,k+£-l 

L  Cr 
pr  =  ~2 


pjj-H,k  "  2pjj,k  +  pjj-l,k 


pjj+l,k  +  2pjj,k  +  pjj-l,k 


where  L  =  n  or  n+1,  t  =  0  or  1-2*NFLIP,  and  jj  =  j  or  j+1, 

depending  on  whether  the  damping  terms  apply  for  the  predictor 

step,  eq.  (17a) ,  or  the  corrector  step,  eq.  (17b) .  These 

damping  terms  are  similar  to  those  used  by  Baldwin  and 

MacCormack  (ref.  11) .  The  pressure  terms  that  they  originally 

used  have  here  been  replaced  with  the  densities .  It  was  found 

that  the  tangential  discontinuity  emanating  from  the  leading 

edge  of  the  wing  forced  large  oscillations  in  the  densities 

and  very  little  in  the  pressures.  Thus  the  original  damping 

did  not  adequately  control  the  nonlinear  instability  near  the 

wing  leading  edges.  The  arbitrary  coefficients  Cr  and  Cm  range 

from  0  for  no  damping  to  1  for  the  maximum  allowable  damping. 

Typical  values  were  <3  =  0.1  and  C  =  0.5,  for  the  radial  and 

r  m 

meridional  damping,  respectively. 

The  above  finite-difference  equations  are  applied  only  to 
the  field  points.  At  the  body  (j  =  3,  3  <.  k  NPHI)  Abbett's 
scheme  is  used  to  satisfy  surface  tangency  as  discussed  later. 
It  requires  information  provided  by  the  finite  difference 
algorithm.  The  algorithm  used  for  the  field  points,  however, 
cannot  be  used  on  the  surface  since  it  requires  points  on  both 
sides  of  the  point  being  advanced  and  thus  data  at  a  set  of 
points  that  would  lie  inside  the  body.  Therefore,  a  second- 
order  accurate  algorithm  is  used  requiring  data  only  on  or 


outside  the  body.  This  scheme  uses  the  predictor  step  of 
MacCormack’s  method,  eq.  (17a),  without  the  fourth-order  damping 
terms,  followed  by  the  corrector  step  given  by 


,n+l 

'‘jr* 


1 

2 

A? 


E1?  .  + 

1/k  },k 


A ,n+l  I  , 1 

45 _  Fn+1 


At 


_  Fn+l' 

j+l,k  j,k 


n+1 


An 

n+l 


,n+l 


rn+l  _  r 
[  j,KCL  j,KCRj 


-  A^n+1 


A? 


At 


(Fj+2,k  2Fj+l,k  +  Fj,k 


(19) 


where  j  =  3.  After  MacCormack's  predictor  and  eq.  (19)  have 
been  used  to  advance  the  data  at  the  body,  Abbett's  scheme  is 
used  as  a  final  corrector.  This  scheme  is  described  in  a  later 
section. 


At  the  shock  wave,  a  predictor-corrector  sequence  is 
again  used  and  requires  data  at  the  shock  and  one  point  below 
it.  The  algorithm  is  as  follows: 


Predictor: 


n+i 

Ej.k 


_  »  _.n+ 1 

e1?  -  — 

3  ,k  At 


„n  _n 
lFj,k  "  j-lfkj 


_  44 


n+l 


An 


rn  -  r 
[j,KPL  j , KPRj 


'  A'n+1  Hj,k 


Corrector: 


n+! 

^j,k 


1 

2 

A? 


E1?  ,  +  E  *}+} 

D»k  i,k 


AC 


n+l 


At 


n+l 


An 


n+l  rn+l 
[  j,KCL  ~  j , KCRj 


„n+l  n+l 
[Fj,k  Fj-l,k 


-  Acn+1 


(20a) 


(20b) 


where  j  =  NT2 .  Equations  (20a)  and  (20b)  are  used  in  conjunc¬ 
tion  with  the  Rankine-Hugoniot  relations,  described  later,  to 
determine  the  peripheral  shock  shape.  It  should  be  noted  that 
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the  damping  terms  are  not  used  for  the  boundary  finite-differ¬ 
ence  equations. 

The  integration  stepsize  is  determined  by  the  stability 
bound  of  the  numerical  scheme.  The  stepsize  must  be  small 
enough  to  satisfy  the  stability  requirement  of  the  scheme  but 
large  enough  to  finish  the  computation  with  a  minimum  of  time. 
The  use  of  the  largest  possible  stepsize  for  hyperbolic  equa¬ 
tions  insures  that  the  finite-difference  scheme  is  as  nearly 
compatible  with  the  method  of  characteristics  as  possible. 

The  bound  on  the  stepsize  value  is  obtained  by  the  amplifica¬ 
tion  matrix  theory.  This  method  is  based  on  a  locally  linear 
(frozen  coefficient)  analysis  of  the  governing  partial  differ¬ 
ential  equations  coupled  with  a  discrete  harmonic  analysis  of 
the  linear  difference  scheme.  The  partial  differential  equa¬ 
tion  system,  eq.  (9) ,  is  rewritten  in  primitive  variable, 
nonconservative  form  as 


(21) 


where  U  =  (u, v,w,p, p) 1 ,  and  M  and  N  are  the  coefficient  matrices 
given  below.  The  amplification  matrix  theory  for  the  two- 
dimensional  space  requires  that 


and 


4^  <  l/(oM) 

At  M  max 


(oM)  =  !  a  (M)  | 


local  max 


(22) 


where  a.,  is  the  local  maximum  modulus  of  the  eigenvalues  of  the 
matrix  M  for  a  given  grid  point  in  the  field.  The  maximum  of 
the  local  maximum  of  the  eigenvalues  in  a  particular  c  =  con¬ 
stant  plane  is  what  is  needed  in  the  above  equation.  A  similar 
condition  is  obtained  in  the  z,n  space 


An 


1/,(cVmax; 


oN  =  |o(N) 


local  max 


(23) 
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This  planar  analysis  has  been  shown  to  give  a  good  bound  on  the 
stepsize  in  multidimensional  problems  if  eqs.  (22)  and  (23)  are 
replaced  by 


—  =  CFL/(oM) 

At  '  M  max 


(24a) 


=  CFL/ (cr  ) 

An  N  max 


(24b) 


where  CFL  is  the  Courant-Friedricks-Lewy  number.  The  CFL 
number  <  1  can  be  varied  during  the  computation  and  is  usually 
assigned  a  value  of  0.9  or  less.  When  applying  eq.  (24)  to 
determine  As,  the  mesh  spacing  Ax  and  An  are  given.  It  is 
therefore  necessary  to  determine  the  minimum  As  predicted  by 
the  two  relations.  This  minimum  As,  which  is  recalculated 
after  every  few  (typically  five  or  less)  integration  steps,  is 
the  one  used  for  the  succeeding  steps  in  the  integration  pro¬ 
cedure. 


The  values  0M  and  cN  required  in  eq.  (24)  are  determined 
from  the  matrices  M  and  N.  The  coefficient  matrix  M  is  given 
by 


M 


ll  i  +  il 

3z  3r 


<a*1b>  +  U  <»‘lo 


(25a) 


and  N  is  given  by 

N  =  h  (A’lc) 

where  I  is  the  identity  matrix  and  (A  ^B)  is  given  by 


(25b) 


where  a  is  the  speed  of  sound  and  is  equal  to  /Yp/ (kp)  in  terms 
of  the  nondimensional  variables  used  in  this  report.  The  five 
eigenvalues  of  M  are 


M  5t  U  1  a/(U2-a2)(m 2+m2)  +  (0-^oT 

,'1,2“S{3I+  02_a2 


(27a) 


M  3t  9£  Q 

3,4,5  3?  9z  U 


(27b) 


where 


The  five  eigenvalues  of  N  are 


|i  V  H 

-  w 

3r 

3<j> 

r 

+  i 

3r  z 

3<J> 

r 

+  , 

3r  y 

r 

+  li , 

9r  x 

9  <f> 

are 

N  _  9n  ^-V2  1  a/(U2-a2)  (<}>2+4>2)  +  (W-^u) 

°1,2  =  3*  ^2 


(28a) 


N  =  3n  .  W 
3,4,5  9<(>  U 


(28b) 


Equations  (27a)  and  (28a)  are  the  slopes  of  the  charac¬ 
teristics  in  the  c,x  and  c,n  planes,  respectively,  while  eqs. 
(27b)  and  (28b)  are  the  slopes  of  the  streamlines  in  their 
respective  planes.  The  maximum  modulus  of  eigenvalues  used 


in  eq.  (24)  is  then  obtained  from  eqs.  (27a)  and  (28a) 


0  ^  1  local  max 

=  max  ( |  a  ^ 1 , 

|o?|) 

(29a 

a (N)  I .  . 

1  local  max 

i  N  i 

=  max(  |  , 

(29b 

5.  BOUNDARY  AND  INITIAL  CONDITIONS 

Several  types  of  boundary  conditions  must  be  satisfied  for 
the  wing-body  problem  considered  in  this  report.  These  include 
solid  surfaces  such  as  the  wing-body  surface,  permeable  sur¬ 
faces  such  as  the  fitted  outer  bow  shock  wave,  symmetry  planes, 
and  initial  planes  from  which  the  computation  can  be  started. 
Another  type  of  boundary  to  consider  is  at  points  on  the  body 
where  the  transformation  becomes  singular;  points  such  as  the 
wing  tip  or  wing-body  juncture.  In  both  of  these  locations  the 
transformation  metrics  vanish  or  become  infinite  and  the  numer¬ 
ical  procedure  breaks  down.  The  separate  procedure  devised  to 
obtain  the  flow  variables  at  these  points  is  described  below. 

A.  Body  Surface 

The  success  of  the  numerical  scheme  depends  in  a  large 
part  on  the  proper  treatment  of  the  surface  boundary  condition. 
For  the  inviscid  flows  considered  in  this  report,  the  fluid 
velocity  vectors  on  the  body  lie  in  the  surface  tangent  planes. 
The  scheme  devised  by  Abbett,  reference  7,  is  used  to  satisfy 
the  surface  tangency  condition.  The  scheme  which  combines 
simplicity  with  accuracy  is  only  applicable  for  steady  super¬ 
sonic  flows.  It  treats  the  computation  at  the  body  in  the  same 
predictor-corrector  fashion  as  the  interior  points  and  is, 
therefore,  compatible  with  the  rest  of  the  calculations. 

Abbett' s  scheme  is  a  predictor-corrector  procedure  in 
which  the  values  of  the  flow  variables  (u,v,w,p,p)  are  first 


predicted  at  =  ?n  +  A 5  using  the  MacCormack  scheme  and 

then  these  values  are  corrected  using  simple  expansion  or  com¬ 
pression  waves  to  enforce  the  surface  tangency  condition 
exactly.  By  restricting  the  flows  to  be  steady  and  supersonic, 
Abbett  is  able  to  use  the  analytical  solution  of  Prandtl-Meyer 
expansion. 


„n+l 


.n 


The  flow  variables  at  i**’’1'  =  5“  +  A£  are  predicted  numer¬ 
ically  by  using  eq.  (17a)  without  the  additional  dissipation 
terms  for  the  predictor  step  and  eq.  (19)  for  the  corrector 
step.  Decoding  the  conservative  variables,  call  them  E^?  from 

(1)  "(1)  -ft.  pft.  - ’ft  ^ 


eq.  (19)  yields  u^,  v3,k'  w’  +■'  p- 


The  body  surface  is  given  by 


f b ( 2 , r , )  =  r  -  rb(z,<{>)  =  0 


-  r  -  1 


(30) 


=  0 


The  unit  outward  normal  vector  is  in  the  left-handed  Cartesian 
coordinate  system 


nb  = 


grad  ffo 
grad  ffa 


-ri-rj-rk 
x  y  x 


(r2  +  r2  +  r2) 
x  y  z 


1/2 


(31) 


The  surface  tangency  condition  to  be  satisfied  on  the  body  is 

q  •  nb  =  0 


or 


Urx  +  Vry  +  WrZ  =  0 


(32) 


with  the  velocity  vector  given  by  its  Cartesian  components. 
The  velocity  vector  at  the  surface 


U3,kX 


+  v 


+  W3,kk 


(33) 


predicted  by  the  numerical  scheme  eqs.  (17a)  and  (19)  will  not 
satisfy  the  tangency  condition  and  will  be  rotated  out  of  the 


tangency  plane  by 


A9  =  sin 


■1  -Ml)  -»•  /  (1) 

_93,k  '  V«3,k_ 


(34) 


where  is  the  modulus  of  and 


-u(1)r  -  v(1)r  -  w(1)r 

.  S  =  U3 ,  krx _ ?j,kfx _ !%jcfz 

q3,k  b 


(35) 


/TT  F~  2 

/r  +  r  +  r 
x  y  z 


and  the  subscripts  3,k  refer  to  the  body  surface.  The  flow  can 
now  be  rotated  according  to  the  Prandtl-Meyer  expansion  or 
compression  given  the  A8.  If  A0  is  positive,  an  expansion  is 
necessary  to  rotate  q^^  into  the  surface  tangent  plane  and 
if  A6  is  negative  a  compression  is  necessary. 

If  the  rotation  angle  A0  is  small  (under  5°),  the  following 
asymptotic  expansion  of  the  Prandtl-Meyer  expansion,  reference 
10,  gives  the  final  corrected  pressure  at  the  surface 


>n+1  =  p(D 
p3,k  p3,k 


Y(M^£)2 

1 - idS -  A0 


/(M 


3,  k 


+  ^(M3^>2 


(V+l)  (M^)4  -  4((M^)2-1] 

4[(M^^)2-1]2 


A0‘ 


+  0  (A0 °)  (36) 


where 


M 


a(1) 
(1)  _  q3,k 

3,k  TuT 

a3,k 


and 


i(1)  = 

3,k 


P(1) 
Y-l  p3,k 

2  (1) 

p3,k 


1/2 
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Assuming  the  entropy  to  remain  constant  during  the  rotation 
•  n+  X 

of  the  velocity  vector,  the  density  ^  can  be  determined  from 


rtn+l  _  .  n+l/v.l/Y 
p3,k  *  (p3,k/K) 


(37) 


with  K  « 

The  velocity  magnitude  is  obtained  from  the  energy  equa¬ 
tion,  eq.  (2) , 

n+1  n+1,  n+1. 1/2 

q3>k  =  U  "  P,  u/P,  ,.) 


3 ,  k  3  ,k 


It  remains  to  find  the  individual  Cartesian  velocity  com¬ 
ponents  of  q"+k-  Since  the  predicted  velocity  vector  was 
rotated  in  the  plane  of  itself  and  the  surface  normal  vector, 
the  direction  of  the  final  velocity  vector  is  n^,  that  is, 


where 


-►n+1  n+1  •+■ 

cr_  .  *  n  •  n 


l3,k 


‘3,k 


n. 


= 


q(1) 

q3,k 


-  (q 


|q(1) 

Iq3,k 


(1) 

3,k 


nb,nb 


-  (q 


(T) 

3,k 


nb)nbl 


The  components  are  now  easily  obtained  as 


(38) 


M 


n+1  _  n+1  _ 

3,k  q3,k  [m2  +  N2  +  L2]1^2 


n+1  _  _n+l  N 

V3,k  -  q3,k  [S2  ;  ?  +  Z2}l/2 


(39) 


n+1  _  n+1  L 

w3,k  -  «3,k  [R2  +  5i  ;  riji/2 
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where 


M  ■  u(1)  - 


N  =  v 


L  =  w 


3,k 

(1) 
3,  k 

(1) 

3,k 


(u3.k'rx  +  v3^k‘ry  *  w3tk-r2)rx 

I  I2 

(u31k’rx  +  V31k'ry  +  w31k*rz)rv 

I  I2 

(u'^.r  +  v'^.r  +  W^’.r  )r 

3 ,  k  x _ 3 ,  k  y _ 3,k  z  z 

i  1 2 


and 


2 

r 


x 


+  r 


+  r 


2 

z 


This  completes  the  boundary  procedure  for  the  body  surface. 
It  is  essentially  similar  to  the  procedure  described  by  Kutler, 
et  al.,  reference  3.  The  most  important  difference,  other  than 
the  different  geometries  considered  and  the  generalized  coor¬ 
dinate  system,  is  the  means  by  which  the  constant,  K,  eg.  (37), 
is  computed.  In  Kutler,  et  al.'s  work  K  is  assumed  to  be 
constant.  For  the  configurations  considered  by  them  a  single 
stagnation  streamline  wetted  the  entire  body,  hence  the  entropy 
on  the  body  is  fixed.  However  the  assumption  of  constant 
entropy  at  the  body  is  violated  if  crossflow  shocks  appear 
which  intersect  the  body.  They  argue  however  that  these  cross- 
flow  shocks  are  weak  and  will  not  significantly  effect  the 
solutions.  Thus  once  the  stagnation  streamline  is  located, 
the  body  surface  entropy  and  thus  K  is  known  immediately  from 
the  flow  variables  on  any  point  on  that  streamline.  Their  work 
is  restricted  to  nonseparated  crossflows.  In  the  present  work 
the  crossflow  must  separate,  and,  hence,  the  body  is  no  longer 
wetted  by  a  single  stagnation  streamline.  The  "constant",  K, 
is  no  longer  constant  over  the  entire  body  since  the  entropy 
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srwpp 


is  not  constant.  K  must  thus  be  determined  numerically  at 
each  point  on  the  body  surface. 

B.  Shock  Surface 

The  boundary  condition  at  the  shock  surface  is  similarly 
treated  as  at  the  body  surface  except  that  the  Rankine-Hugoniot 
conditions  must  be  satisfied  instead  of  the  flow  tangency  con¬ 
dition.  Specifically  we  need  to  determine  the  shock  surface. 
Let  the  shock  surface  be  given  by 


f  (z,r,<t>)  =  r-r  (z ,  <J>)  =  0 

S  o 


(40) 


The  outward  unit  normal  is  given  in  terms  of  the  left-handed 
Cartesian  unit  vectors  by  (see  fig.  6) 


nx 


yyy  {(r*'rs**x)i +  <vrs*Vj+  “v 


rsrrs4,*z»k 


(41) 


where 


Vf 


'  1  (rx-rs4.*xl 


+  Vrs*V2  +  lvI.rr.t,.|!| 


1/2 


The  metrics  r  ,r  ,r  ,  <J>  ,$  ,  and  <J>  are  known  once  r  and  <J>  of 

A  V  Z  A  J  “ 

the  shock  surface  are  known.  The  preshock  velocity  vectors 
normal  and  tangential  to  the  shock  surface  are 


and 


us,ns 


(42) 


v1  =  q  -  u.  ,  respectively 

X  00  X 


(43) 


The  notation  used  in  this  section  is  as  follows:  Subscript  1 
or  00  denotes  preshock  condition,  subscript  2  -  postshock, 
superscript  tildas  on  u  and  v  denote  the  velocity  vectors 
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normal  and  tangential  to  the  shock  surface,  respectively.  The 
magnitude  of  the  normal  preshock  velocity  vector  is 


u-,  =  (q  *n_)  /  n_*n. 


100  s  s  s 


Vf 


b  (,r-‘ 


r  ,  ) u  +  (r-r  .<t>  )v 

s<)r  x  00  y  s<}>yy'  ° 


+  (r  -r  -r  .  4>  ) 
z  sz  s<tyz 


Vf  I  (44) 
s  1 


Squaring  and  cancelling  obtains 


(r  -r  ,<p  )  u  +  (r  -r  ,4  )v  +  (r  -r  -r  <6  ) 
x  si  x  °°  y  sifry'  °°  v  z  sz  s<}rz; 


w“}: 


Solve  now  for  r  to  obtain 
sz 


rsz  =  (rz“rs^z)  +  w»n  +  51^2  +  { 


(45) 


where 


„  .  |rX  -  *  lrV  ~  r5tl>V>V. 


2  ~  2 
w  -  u. 

00  1 


(  >  ‘  (rx  -  rs**x>2  +  (ry  -  rs*V2 


The  shock  wave  is  moved  according  to  the  following  Euler 
predictor  and  modified  Euler  corrector: 


Predictor: 

rn+1  = 
s 

n  . 
r  + 
s 

n 

rsz 

4;n+1 

(46) 

Corrector : 

rn+1  = 
s 

H 

rn  +  r 
sz 

(n+1) 

sz 

4t"+1 

(47) 

34 


where  A£  is  the  integration  stepsize.  The  quantities  rg^  are 
evaluated  by  simple  one-sided  differencing 


r 


n+1 

SNT2 / KPL 


n+1 


NT2, KPR 


/ 


PNT2 , KPL 


NT 2, KPR 


(485 


and 


r 


n+1 

S<}> 


n+1 

SNT2 , KCL 


s+1  ]/ 

NT2,KCR' 


*NT2,KCL 


^NT2 , KCRj 


(49) 


The  subscript  indices  have  been  previously  defined.  The  pro¬ 
cedure  for  determining  the  shock  slope  is  as  follows:  the 
postshock  conditions  are  predicted  using  eq.  (20a)  and  eq.  (46) . 
Decoding  obtains  a  predicted  postshock  pressure,  p2,  from  which 
the  normal  preshock  velocity  is  given  for  a  perfect  gas  by 


Pa 


Pj 


p°°  y2-1  p2 

2pco  L  2Y  Pco 


(Y-l) 


(50) 


Forming  the  metrics  from  the  predicted  rg+^  and  rj^1  from 
eq.  (48)  gives  sufficient  information  to  solve  for  the  shock 
slope  using  eq.  (45) .  The  postshock  velocity  vector  can  now 
be  determined  from 


->• 


(51) 


where  u2  and  v2  are  the  postshock  normal  and  tangential  velocity 
components.  The  tangential  component  remains  unchanged  across 
the  shock,  thus 


(52) 


and  the  normal  component  is  determined  by  use  of  the  con¬ 
tinuity  equation 

S2  *  Gipi/p2  (53) 
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and  p2  is  given  by  the  Rankine-Hugoniot  condition 


p2(p2;  Pco'P~>  =  P 


oo  00 


£2  Y-l 
Poo  Y+l 

1  +  Izi  £2 

TT+1  P_ 


(54) 


The  normal  postshock  velocity  is  given  by 


q,  =  u_i  +  v0j  +  w_k 


=  qn 


-  fi .  fii 

l  »,) 


u,  n 
1  s 


(55) 


which  can  be  broken  into  the  requisite  Cartesian  postshock 
velocity  components  u2,v2,w2.  This  completes  the  predictor 
step.  The  entire  procedure  is  repeated  for  the  corrector  step 
except  that  the  appropriate  equations  to  be  used  are  now 
(20b),  (47),  and  (49). 


C.  Symmetry  Plane 

The  flow  fields  of  interest  in  this  report  have  a  plane 
of  bilateral  symmetry  since  no  yaw  angles  are  considered  and 
the  configurations  considered  are  bilaterally  symmetric.  This 
allowed  us  to  reduce  the  number  of  mesh  points  required.  The 
symmetry  conditions  are  imposed  by  simple  reflection  of  the 
flow  variables  and  the  outer  bow  shock  location  about  the  two 
planes  of  symmetry  4>  =  0°  and  <p  =  180°.  Since  MacCormack’s 
scheme  is  a  noncentered  predictor-corrector  scheme,  it  is  not 
possible  to  apply  the  symmetry  conditions  after  the  predictor 
step.  The  predicted  values  of  the  flow  variables  are  obtained 
directly  by  applying  the  predictor  step  eq.  (17a) ,  to  the  point 
beyond  the  symmetry  planes,  i.e.  <p  =  0°  -  A<f>  and  <J>  =  180°  +  A$. 
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The  symmetry  conditions  are  imposed  only  after  the  corrector 
step  to  provide  sufficient  information  for  the  predictor  step 
for  the  next  marching  step. 

Treating  the  symmetry  planes  in  this  manner  avoids  the 
spurious  oscillations  emanating  from  the  symmetry  planes. 
These  oscillations  occur  when  MacCormack's  scheme  is  employed 
with  a  generalized  curvilinear  coordinate  system  and  symmetry 
is  imposed  after  both  the  predictor  and  corrector  steps. 


D.  Wing-Tip  Condition 

At  the  wing-tip  the  transformed  Euler  equations,  eq.  (6) 
are  indeterminate.  This  can  be  snown  as  follows:  the  basic 
equations  are  eqs.  (4)  along  with  the  contravariant  velocity 
components,  eqs.  (5) .  The  appropriate  metrics  are  given  by 
eqs.  (6) .  If  we  now  consider  the  region  near  the  wing  tip 

r  =  1  +  e  J  e |  <<  1 

and  (56) 

e  ~  e  | e ]  <<  l 


where  4>  -  ir/ 2  =  0 ,  we  can  write  the  metrics  in  asymptotic  form 
as 


xr  =  ae 


x0  =  -a0 


yr  =  ae 


Y0  =  ae 


where 


1  ± 


A  -  4B2/R2 


and  the  Jacobian  is 


-1-22  2 
J  =  a  (e  +  0  j 


(57) 


(58) 


(59) 


With  these  metrics,  the  asymptotic  form  of  the  governing 
equations  becomes 

PeV  +  P0W  +  PC  +  pD  =  0 

(pu) £V  +  (pu)QW  +  pu (C+D)  =  0 

(pv) £V  +  (pv)QW  +  pv (C+D)  =  0 

(pw) ev  +  (Pw)ew  +  pw (c+d)  +  e  kp  +  e  kp  =  0 

t.  Z  u 


where  the  overbar  (  )  stands  for  J-1 (  ) 

C  =  a [eu£  +  ev£  +  we(-6y2  +  ex  )  ] 

and 

D  =  a[-0uQ  +  evQ  +  w0(xze  -  yze)] 

All  of  the  overbarred  quantities  vanish  as  e,0  ■>  0.  Thus  as 
the  wing  tip  is  approached,  the  basic  equations  are  always 
satisfied  provided  the  flow  variables  are  continuous.  Thus 
the  Euler  equations  can  be  satisfied  by  any  continuous  values 
of  the  primitive  variables  at  the  wing  tip. 

If  we  now  allow  for  discontinuous  solutions  of  the  Euler 
equations  at  the  wing  tip  we  have 

[fi]  cos(n,z)  +  [F]  cos{n,r)  +  [G]  cos(n,<J))  =  0 

where  the  cos(n,z),  cos(n,r),  and  cos(n,<|>)  are  direction 
cosines  between  the  coordinates  and  the  surface  of  discontin¬ 
uity.  From  the  previous  asymptotic  forms  of  the  metrics,  the 
jump  relations  to  the  lowest  order  become 
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[pV]  cos(n,r)  +  [pW]  cos(n,<j>)  =  0 

[  puV  +  y0kp]  cos  (n,r)  +  [puW  -  yrkp]  cos(n,iji)  =  0 

[pvV  -  xQkp]  cos  (n, r)  +  [pvW  +  xrkp]  cos(n ,  $)  =  0 

[pwW  +  e2kp]  cos(n,r)  +  [pwW  +  6zkp]  cos(n,4>)  =  0 


Since  all  the  metrics  and  the  overbarred  quantities  vanish  at 
the  tip,  it  follows  that  the  jump  relations  are  identically 
satisfied  there. 

Thus  for  both  discontinuous  and  continuous  flows,  the 
Euler  equations  are  indeterminate  at  the  wing  tip.  In  other 
words  arbitrary  values  can  be  set  for  p,u,v,w  at  the  wing  tip. 

Since  we  are  free  to  pick  the  primitive  flow  variables 
at  the  wing  tip,  we  can  be  guided  by  the  actual  physics  of  the 
flow.  We  will  model  the  flow  of  a  vortex  sheet  coming  off 
tangentially  to  the  wing  at  the  tip.  It  is  well  known  that 
for  compressible  flows  that  tangency  cannot  be  satisfied  for 
positive  pressures  and  densities  if  the  flow  must  turn  180° 
as  at  the  wing  tip.  Thus  we  require  the  flow  to  separate. 

Even  though  the  equations  are  indeterminate  at  the  wing  tip, 
the  primitive  variables  must  satisfy  the  jump  relations  for 
the  untransformed  Euler  equations  and  these  are  for  a  tangen¬ 
tial  vortex  sheet 


(i)  p  is  continuous 

(ii)  v  =  0  (i.e.  the  sheet  is  a  streamline) 

(iii)  total  enthalpy  is  constant 


fT~  2 

/u  +  w 


across  the  sheet 


The  jump  in  the  tangential  velocity 
is  arbitrary.  The  jump  in  density  is  then  determined  from  the 
constancy  of  total  enthalpy. 
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Two  more  conditions  are  required  before  all  the  flow 
variables  are  determined  at  the  wing  tip. 

The  two  additional  conditions  that  we  have  used  with 
partial  success  are: 

(iv)  u/w  =  tan(iji  - 

(v)  p/pY  =  constant 

where  ip  is  the  leading-edge  sweep  angle.  In  particular  we  set 
=  'i'jjj/  i.e.  the  flow  at  the  tip  follows  the  mean  flow  direc¬ 
tion  obtained  by  averaging  the  flow  directions  just  inboard  of 
the  tip  on  the  upper  and  lower  surfaces.  These  five  conditions 
are  similar  to  those  utilized  by  Minailos  (ref.  12)  except  that 
he  set  so  as  to  obtain  good  agreement  with  experimental  data. 

These  conditions  are  implemented  into  the  code  as  follows: 

a)  pQ  =  0.5^  +  p_1)  b)  pQ  -  0.5^  +  p_^) 

c)  vq  =  0  d)  \pQ  =  0.5(1^  + 

e)  total  enthalpy  is  constant 

where  subscript  "o"  designates  the  wing  tip,  "1”  the  first  point 
on  the  upper  surface  inboard  of  the  tip  and  "-1"  the  first  point 
on  the  lower  surface  inboard  of  the  wing  tip.  In  other  words 
the  pressure,  density,  and  flow  angle  at  the  wing  tip  are 
obtained  by  averaging  over  the  nearest  surface  neighboring 
values.  The  reason  for  averaging  the  density  instead  of  using 
condition  (v)  is  that  the  latter  condition  requires  an  accurate 
evaluation  of  ,the  surface  entropy.  This  was  not  possible  with 
the  present  code.  An  alternate  boundary  condition  for  deter¬ 
mining  the  densities  on  each  side  of  the  vortex  sheet  would  be 
to  specify  the  rate  of  vortex  shedding  from  the  leading  edge. 
From  this  information  the  velocity  on  each  side  of  the  sheet 
would  be  given  and  the  densities  could  then  be  calculated. 

This  is  unlike  two-dimensional  airfoil  theory  where  conditions 
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(i) ,  (ii)  and  (iii)  fix  the  strength  of  the  vortex  sheet.  The 
reason  for  this  is  that  in  two-dimensional  flows  the  wake  con¬ 
sists  of  a  single  split  streamline  whereas  for  three-dimensions 
the  vortex  sheet  consists  of  separate  streamlines  from  the  top 
and  bottom  surfaces  which  may  be  of  different  entropies. 


E.  Wing- Body  Juncture 

At  the  wing-body  juncture,  the  transformation  metrics 
(eq.  (16))  become  singular  since  the  term  under  the  radical 
sign  vanishes. 

This  singular  behavior  is  due  to  the  abrupt  jump  of  the 
surface  slope  between  the  body  and  the  wing.  As  mentioned  in 
section  3,  this  singular  behavior  can  be  automatically  avoided 
if  the  analytical  expressions,  eq.  (16) ,  for  the  metrics  are 
replaced  by  their  finite  difference  equivalent.  No  special 
treatment  will  then  be  required  at  these  points  as  was  re¬ 
quired  for  the  singularity  at  the  wing- tip. 

However  to  replace  the  analytic  expressions  by  their 
finite  difference  expression  is  impractical  in  context  of  the 
existing  Space  Shuttle  code,  refs.  3  and  4.  In  that  code  the 
transformation  of  the  physical  body  to  a  unit  radius  cylinder 
as  given  by  eq.  (16)  and  the  normalization  of  the  distance 
between  the  body  surface  and  shock  surface  as  given  by  eq.  (7) 
are  treated  as  two  separate  transformations.  The  mesh  point 
coordinates  from  which  the  finite  differences  would  be  com¬ 
puted  are  determined  only  as  the  net  result  of  both  transforma¬ 
tions.  Thus  there  is  no  straightforward  procedure  for  comput¬ 
ing  the  metrics,  eq.  (16),  numerically. 

The  procedure  used  in  this  report  is  to  avoid  the  wing- 
body  juncture  singularity.  This  is  done  by  filleting  the 
juncture  to  eliminate  the  abrupt  change  in  surface  slope.  The 
simplest  method  of  achieving  this  is  to  add  a  small  perturbation 
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[ 

to  the  transformed  body  radius,  so  that  instead  of  r^z,*!))  =  1, 
we  have  r ^  ( z ,  <}> )  =  1  +  e.  The  perturbation,  e,  used  herein  is 
0.025.  With  this  perturbation  the  wing-body  configuration 
shown  in  figure  5a  is  obtained. 

F.  Initial  Starting  Plane  Conditions 

The  starting  solution  or  initial  data  plane  can  be  obtained 
by  one  of  two  ways.  The  Space  Shuttle  code  has  the  capability 
to  generate  internally  the  solution  for  a  pointed  cone  at  angle 
of  attack  using  the  shock-fitting  procedure.  To  start  that 
solution  all  the  grid  points  in  the  initial  plane  are  initial¬ 
ized  using  the  analytical  solution  for  a  pointed  cone  at  zero 
degree  angle-of -attack  and  the  given  free-stream  Mach  number. 

The  integration  is  started  and  the  angle  of  attack  is  incremented 
over  100  marching  steps  to  its  final  value.  Thereafter  the 
integration  is  continued  until  the  calculation  has  converged 
to  a  conical  flow  field  solution.  Since  the  flow  is  conical, 
the  axial  coordinate  is  stepped  back  to  its  original  starting 
value  to  give  the  starting  plane  flow  field  for  a  pointed  cone  at 
angle-of-attack.  This  procedure  was  used  to  provide  the  initial 
data  for  the  wing-body  configurations  studied  in  this  report. 

For  the  delta  wing  alone,  a  slightly  different  procedure 
is  required.  Starting  with  a  given  wing-body  solution  of  the 
desired  sweep  angle,  angle-of-attack,  and  free-stream  Mach 
number,  the  integration  procedure  is  continued  while  the 
physical  body  radius  is  gradually  forced  to  vanish  over  twenty 
to  thirty  marching  steps.  Thereafter  the  transformed  body 
radius,  r^(z,<p),  is  forced  slowly  to  unity  over  another  twenty 
steps.  The  integration  continues  until  again  a  conical  flow 
field  is  obtained.  If  solutions  are  required  for  planar  wings 
other  than  delta  wings,  then  the  conical  planar  delta  wing 
solution  is  used  as  a  starting  solution  and  the  sweep  angle  can 
vary  arbitrarily  with  the  axial  coordinate,  z. 
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6.  NUMERICAL  RESULTS  AND  DISCUSSION 

Numerical  calculations  have  been  obtained  for  two  cases: 
a  planar  delta  wing  alone  of  aspect  ratio  1.0  at  free-stream 
Mach  number  3  and  angle  of  attack  a  =  10° ,  and  a  delta 
wing  in  the  presence  of  a  circular  body  at  Mb  =  3  and  a  =  10°. 

A.  Wing  Alone 

The  first  numerical  results  are  for  delta  wing  alone  at 
M^  =  3  and  a  =  10°.  The  closest  available  experimental  data, 
reference  13,  are  for  a  configuration  shown  in  figure  7(a). 

This  experimental  configuration  has  a  finite  thickness  delta 
wing  with  the  edges  beveled  to  obtain  sharp  leading  and  trail¬ 
ing  edges.  The  experimental  free-stream  Mach  number  is  2.86. 
The  pressure  coefficients  are  shown  in  figure  7 (b) .  The 
effects  of  thickness  have  been  subtracted  out  of  the  experi¬ 
mental  results  shown.  The  present  Euler  code  results  and  the 
experimental  data  compare  quite  well  except  near  the  tip.  The 
effects  of  the  bevel  at  the  tip  of  the  experimental  configura¬ 
tion  cause  some  scatter  in  the  data.  The  cause  of  the  spikes 
in  the  numerical  pressure  coefficients  is  not  clear.  It  is 
suspected  that  they  are  a  result  of  the  nonlinear  instability 
of  MacCormack's  scheme  at  the  tangential  discontinuity  at  the 
wing  tip.  The  lift  coefficients  are  0.222  and  0.215  for  the 
Euler  code  and  experimental  results,  respectively. 

The  conical  streamline  plot,  figure  8,  shows  the  location 
of  the  two  vortical  singularities  and  the  stagnation  stream 
surface.  Vortical  singularities  occur  at  the  windward  plane 
of  symmetry  and  above  the  upper  surface  30%  of  the  span  inboard 
of  the  leading  edge.  The  stagnation  stream  surface  is  on  the 
lower  surface  about  1%  span  inboard  of  the  wing  tip.  These 
results  show  that  the  Kutta  condition  imposed  at  the  wing  lead¬ 
ing  edge  cause  the  flow  to  separate  from  the  tip  tangentially 
to  the  surface.  There  are  no  experimental  data  such  as  vapor 
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screen  data  for  this  angle  of  attack  and  Mach  number  to  verify 
the  location  of  the  leading-edge  vortex  for  a  thin  delta  wing. 
Figure  9  shows  the  particle  paths  obtained  by  the  numerical 
results  as  would  be  seen  by  a  vapor-screen  visualization.  The 
flow  field  for  the  delta  wing  are  shown  in  figure  10.  Shown 
are  the  Cartesian  crossflow  velocity  vectors. 

B.  Wing-Body 

The  results  for  the  wing-body  calculation  are  shown  in 
figures  11(a)  to  (c) ,  in  terms  of  pressure  coefficients  at 
stations  just  before  the  wing,  at  the  wing  trailing  edge,  and 
midway  between  these  two.  Figure  11(a)  shows  the  rapid 
expansion  around  the  windward  side  of  the  cylinder  terminating 
with  a  strong  crossflow  shock  wave  followed  by  some  isentropic 
compression  on  the  leeside.  Figure  11(b)  shows  that  the  cross- 
flow  initially  expands  around  the  body  but  then  must  recompress 
to  meet  the  wing-body  juncture.  The  flow  continues  to  compress 
as  the  crossflow  stagnation  point  is  reached  on  the  windward 
side  of  the  wing  just  inboard  of  the  leading  edge.  Rapid 
expansion  around  the  leading  edge  follows  and  then  gradual 
recompression  on  the  leeside  of  the  wing  and  body.  Figure  11(c) 
showing  the  pressure  distribution  at  the  trailing  edge  of  the 
wing  shows  a  continuation  of  these  events.  The  lift  coefficient 
based  on  the  wing  planform  area  of  the  wing  panel  alone  in  the 
presence  of  the  body  is  0.237. 

The  large  oscillations  of  the  pressure  coefficient  seen 
near  the  juncture  of  the  wing  and  body  seem  disturbing  at  first. 
However  their  origin  can  be  pinpointed  to  the  poor  mesh  resolu¬ 
tion  at  the  juncture.  This  is  a  basic  fault  of  any  type  of 
conformal  transformations  such  as  those  used  in  this  report. 
Coordinate  lines  are  quite  strongly  concentrated  at  convex 
corners  such  as  the  wing  tip  and  quite  dispersed  at  concave 
corners.  This  can  be  readily  observed  in  figure  5.  This 
problem  can  be  easily  solved  by  obtaining  the  mesh  points  by 
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some  other  means  than  the  conformal  transformation.  However 
this  will  require  writing  essentially  a  new  computer  code  so 
that  the  transformation  metrics  may  be  computed  numerically. 

As  mentioned  previously  this  is  incompatible  with  the  present 
code. 

The  particle  paths  as  would  be  seen  by  vapor  screen  visu¬ 
alization  is  shown  in  figure  12.  Again  the  flow  is  seen  to  be 
leaving  the  leading  edge  tangentially  to  the  lower  surface. 

C.  Wing-Body  Interference 

The  final  figure  shows  the  wing-body  interference  factor, 
1^,  obtained  experimentally  (ref.  1)  and  by  the  present  results. 
This  factor  is  the  ratio  of  the  lift  on  the  wing  panels  mounted 
on  the  body  to  the  lift  of  the  wing  alone  for  the  same  angles 
of  attack.  Slender-body  theory  indicates  that  a  favorable 
interference  factor  of  =  1.45  is  obtained  on  the  wing  lift 
by  the  presence  of  the  body.  Compressibility  effects  reduce 
this  favorable  interference  to  1.07  obtained  by  the  present 
Euler  code  and  Kutta  condition  compared  to  0.98  experimentally. 
Although  the  airfoil  section  differ  for  the  numerical  and 
experimental  results,  the  adverse  effect  of  compressibility 
on  wing-body  interference  is  predicted  fairly  well. 


7.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  specific  results  obtained  in  this  study  are  as 
follows: 

a)  The  lifting  coefficients  and  flow  field  data  of  a 
M^  =  3  planar  delta  wing  at  an  angle-of-attack 
a  =  10°  have  been  obtained.  The  normal-force  coef¬ 
ficient  compares  well  with  the  experimental 

coefficient  (at  M  =2.86  and  finite  thickness  wing) 

00 

of  0.215. 
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b)  The  normal-force  coefficient,  flow  field  data,  and  the 
interference  coefficient  have  been  obtained  for  a 
delta  wing  in  the  presence  of  a  circular  body  at 
Mm  =  3.0  and  a  =  10°.  The  normal-force  coefficient 
based  on  the  wing  panel  planform  area  is  0.237.  The 
interference  coefficient  Kw  =  1.07  as  compared  to  the 
experimental  Kw  =  0.98  for  a  wing-body  configuration 
with  slightly  different  airfoil  sections  for  the  wing 
panel. 

A  number  of  conclusions  may  be  drawn  from  this  study. 

Good  flow  field  and  surface  pressure  results  are  obtainable 
for  wings  alone  and  wing-body  configurations  with  sharp  leading 
edges  by  combining  an  Euler  equation  solver  with  a  Kutta  con¬ 
dition  to  simulate  the  viscous  effects  at  the  leading  edges. 

The  separation  and  leading-edge  vortex  seem  to  be  properly 
predicted.  However  no  flow  field  comparisons  with  experimental 
data  were  made  due  to  the  lack  of  experimental  flow  field  data. 
Presently  the  speed  of  the  code  is  controlled  by  the  fine  mesh 
at  the  tip  and  the  accuracy  by  the  coarse  mesh  at  the  juncture. 
Thus  the  numerical  procedure  can  be  considerably  speeded  up  and 
the  accuracy  substantially  improved  by  a  better  distribution  of 
mesh  points.  The  additional  dissipation  appended  to  MacCormack's 
scheme  to  control  the  severe  nonlinear  instability  near  the  wing 
leading  edge  is  sufficient  to  stabilize  the  scheme.  But  it  is 
not  known  whether  it  may  be  too  dissipative  and  diffuse  the 
vorticity  too  rapidly.  It  should  be  noted  that  the  application 
of  the  Kutta  condition  is  also  possible  for  flows  with  subsonic 
axial  Mach  numbers. 

Two  recommendations  seem  appropriate.  The  code  should  be 
revised  so  that  the  metrics  are  obtained  numerically  to  allow 
for  a  more  arbitrary  clustering  of  mesh  points  where  necessary, 
i.e.,  more  at  the  wing-body  juncture  and  fewer  at  the  wing  tip. 
The  conformal  transformation  presently  used  clusters  too  much 
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at  the  tip  and  disperses  too  much  at  the  juncture.  Also  a  more 
careful  and  thorough  study  is  required  to  determine  the  best 
stabilizing  terms  to  be  added  to  MacCormack's  scheme.  This 
dissipation  should  be  localized  to  the  region  where  needed. 
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APPENDIX  A 


CARTESIAN  AND  CONICAL  PARTICLE  PATHS 


A  detailed  discussion  is  given  in  this  appendix  on  the 
method  of  obtaining  the  Cartesian  and  conical  particle  paths 
projected  on  a  constant  2-plane. 


Cartesian  Streamline/Particle  Path  Plots 

The  Cartesian  crossflow  velocity  vectors  projected  on  the 
constant  z-plane  trace  out  "pseudo"  streamlines.  These  stream¬ 
lines  represent  steady  state  path  lines  of  fluid  particles. 

The  path  lines  are  determined  by  solving  the  differential 
equation 

^jr  =  u  (x,y)  (A-l) 

af  =  v (x,y)  (A-2) 


A  modified  Euler-Cauchy  method  is  used  to  solve  these 
equations  in  the  form 

t+At 

Ax  =  x*  -  Xj  ^  =  |  u(x,y)dt  (A-3) 

t 

t+At 

Ay  =  y*  -  y_j  k  =  j  v(x,y)dt  (A-4) 

t 


where  x*  and  y*  are  the  projected  locations  of  the  next  point 
on  the  streamline  curve.  The  integration  is  started  after  an 
incremental  time  step,  At,  has  been  defined  and  the  initial 
position  of  the  fluid  particles  at  arbitrarily  selected  point 
have  been  established. 
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The  numerical  integration  is  started  by  initially  assuming 
that  u  and  v  at  (j,k)  are  constant  over  the  interval 

xj  ,k  *  x  *  xj+l,k+l  and  *j,k  *  y  *  yj+l,k+l-  This  results  in 
a  preliminary  estimate  of  x*  and  y*  of 


*  =  xj,k  +  UAt 

(A-5) 

*  *  yj,k  +  v4t 

(A-6 ) 

The  velocity  components  at  (j*,k*)  are  determined  by  linear 
interpolation  as  illustrated  in  the  sketch  below.  We  need  to 


j+1,  k+1 


determine  Ax  and  A^  which  are  defined  by  any  of  eqs.  (A-9) 
and  (A-10) ,  respectively.  We  have 


(A- 7 ) 


(A- 8) 


where 


Ax^xj,k  +  Ax  x j+1 , k 
Ax)yj,k  +  Ax  yj+i#k 
Ax^xj,k+l  +  ^x  xj+l,k+l 
Ax^j,k+1  +  Ax  yj+l,k+l 


xc  =  (1  "  Ay)xj,k  +  Ay  xj,k+l 

yc  =  (1  -  \)yj/k  +  Ay  yj#k+i 

Xd  =  (1  Ay)xj+l,k  +  Ay  Xj+l,k+l 

yd  =  (i  -  yyj+i,k  +  Ay  yj+i,k+i 

Solving  eqs.  (A-7)  and  (A-9)  for  A^  obtains 


Ax  " 


~3x  +  /6x  '  4ax^x 


where 


(A-9) 


(A-10) 


(A-ll ) 


<xj+l,k  "  xj,kf  <yj,k+l  "  yj+l,k+l) 
<xj,k+l  *  xj+l,k+l) (yj+l,k  “  yj,k) 


6x  *  -<xj*l,k  -  xj,k> <yj,k+l  '  yj,k> 

<xj«,k*  "  xj,k)  <yj+l,k  *  yj,k  *  yj+l,k  ”  yj+l,k+l) 
+  (yj+l,k  *  yj,k)<xj,k+l  -  xj,k> 

<yj*,k*  "  yj,kHxj+l,k  ‘  xj,k  '  xj+l,k+l + Xj,k+1; 


Tx  -  (xj*,k*  ‘  xj,k)<yj,k+l  '  yj,k> 


*  (yj*,k*  '  yj,k)(xj,k+l  *  Xj,k> 


i 
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Similarly  eqs.  (A-8)  and  (A-10)  are  solved  for  to  obtain 


Ay  " 


-e„  +  /B*  -  4a  Y. 


(A-13) 


where  a  ,  $  ,  and  y  are  the  same  as  a  ,  6  ,  y  with  the  points 
y  y  y  x  x  x 

(j,k+l)  and  (j+l,k)  interchanged. 

The  interpolated  velocity  components  are  now  easily 
obtained  as 


Uj*,k*  =  (1  ‘  A-)(1  “  AJU-4  .V  +  -  AJu4 


-»  V  VA  «  «  ■  LA  \ 

y  j  ,  k  x 


y  j+lfk 


+  AxAy  Uj+l#k+l  +  (1  Ax)Ay  uj,k+l 


(A- 14) 


vj*,k*  *  (1  '  Ax)(1  -  VVjA  +  V1  ‘  Ay)vj+l,k 
+  Vy  Vl.k+1  +  (1  '  Ax)Ay  Vj,k-H 


(A-15) 


By  averaging  the  appropriate  velocity  components  at 
(j*,k*)  with  those  at  (j,k),  an  improved  estimate  of  x*  and 
y*  is  obtained.  This  improved  estimate  is 


j.k  +  i(uj,k  +  “j..k*)At 

(A-16 ) 

j,k  +  !(vj,k  +  vj*,k*)At 

(A-17) 

The  velocities  at  the  new  (j*,k*)  are  recomputed,  as  above, 
and  then  the  new  location  and  velocity  vector  are  used  as  the 
initial  conditions  for  calculating  the  next  point.  The  process 
is  repeated  until  the  maximum  number  of  time  steps  is  reached. 


Conical  Streamline/Particle  Path  Plot 


For  conical  bodies  such  as  delta  wings  it  is  more  conven¬ 
ient  to  use  conical  streamlines  than  the  Cartesian  crossflow 
streamlines  as  developed  above.  It  is  a  simple  matter  to 
modify  the  above  procedure  to  obtain  the  conical  streamline 
plot.  The  procedure  detailed  below  will  not  give  a  strictly 
correct  streamline  crossflow  but  the  error  is  negligible. 

Suppose  the  vertex  of  the  conical  body  is  at  Zq  and  that 
the  crossflow  data  are  available  at  Zp  as  shown  in  the  sketch 
below. 


At  each  point  (x,y)  in  the  Zp-plane  the  three  Cartesian 
velocity  components  are  given.  We  do  not  have  the  data  given 
on  a  constant  R  surface  as  required  for  a  true  conical  projec¬ 
tion.  However  the  errors  will  be  small  for  -  Z_  >>  r.  The 

P  ° 

required  relations  for  r,  0,  and  4>  are  given  by 

r  =  /x  ^  +  y2 

e  =  tan-1  [r/  (Zp  -  ZQ)  ]  (A-18) 

<j>  =  tan”1  [x/y] 
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We  require  the  (x,y)  components  of  the  conical  crossflow 

velocity  components  so  that  the  procedure  of  the  previous 

section  can  be  employed.  These  components  are  derived  by 

resolving  the  w  component  into  a  component,  call  it  u  , 

c 

along  the  conical  ray  R  and  a  component,  q  ,  in  the  Z  -plane 

c  P 

as  shown  in  the  sketch  below. 


These  components  are  given  by 


u  =  w/cos  9 

G 


q„  =  w  tan  0 
c 


The  (x,y)  components  of  q  are  qiven  by  (see  sketch  below) 

C 

qcx  =  *  %  sin  ♦  4  f 


qcy  =  -  qc  cos  4> 


Adding  these  components  to  the  Cartesian  crossflow  veloc¬ 
ity  components  obtains  the  Cartesian  components  of  the  conical 
crossflow  velocity  vector  projected  on  the  constant  Z  -plane. 


(A-21) 


u  =  u  -  q  sin  4> 
c  c 

v  =  v  -  q  cos  S 
c  c 

These  velocity  components  are  to  be  used  in  place  of  u,v  of 
the  previous  section  to  obtain  the  conical  streamline  path, 
an  example  of  which  is  shown  in  figure  8. 


fRSCIDIAG  PAGS  BUNK-NOT  TILJHED 


i 


B 


C 


a 

a 

B 


C,  D 


CFL 

E,F,G 


(E  ,F,  G,  H) 


LIST  OF  SYMBOLS 


aspect  ratio 

matrix  defined  by  eq.  (26a) 

matrix  defined  by  eq.  (26b) 

speed  of  sound,  Ap/kp 

parameter  defined  by  eq.  (58) 

body  radius  in  Cartesian  space 

parameters  defined  following  eq.  (59) 

damping  coefficients,  defined  following  eq.  (18) 

Courant-Friedricks-Lewy  number 

column  vectors  in  Cartesian  space  (eq.  (1)) 

column  vectors  in  transformed  space  (eq.  (4)) 

column  vectors  in  computational  space  (eq.  (9)) 

components  of  E 

clustering  function  in  meridional  direction  (eq.  (7)  ) 
function  describing  the  transformed  body  surface 
function  describing  the  transformed  shock  surface 
=  x  +  iy 

clustering  function  in  radial  direction  (eq.  (7) ) 

=  AT 

Cartesian  unit  vectors 
identity  matrix 
Jacobian 

indices  in  the  computational  mesh  system  (fig.  5) 
constant  defined  by  eq.  (37) 


( 
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interference  factor 


KCL, KCR 
KPL,KPR 
k 

L,  M,  N 
M 


M 


m  ,  m  ,  m 
x  y  z 

N 

NFLIP 


NPHI 

NT2 


n 

-* 

n 


P 

Q 


q 

R 


u,v,w 


(u,  v,w) 
(x,y,z) 


index  parameters,  defined  by  eq.  (17) 
index  parameters,  defined  by  eq.  (17) 

=  2Y/(Y-1);  or  index  in  computational  space 
parameters  defined  by  eq.  (39) 
coefficient  matrix,  eq.  (25a) 
free-stream  Mach  number 

intermediate  variables  defined  by  eq.  (27) 
coefficient  matrix,  eq.  (25b) 

switching  parameter  defined  following  eq.  (17) 

number  of  meridional  points  (fig.  5) 

number  of  radial  points  (fig.  5) 

index  in  marching  direction 

unit  normal  vector 

pressure 

intermediate  variable  defined  by  eq.  (27) 
velocity,  (/uz  +  vz  +  wz ) 

=  Rw  +  B2/R 

wing  semispan  in  Cartesian  space 

body  radius  in  transformed  space 

shock  radius  in  transformed  space 
T 

=  (u,v,w,p, p) 

contravariant  velocity  components  (eq.  (5)) 

Cartesian  velocity  components  nondimensionalized 
by  maximum  flow  velocity 

Cartesian  coordinates 


z 

(z,r,<p) 

a 

P 

y 

<P 

s 

u 
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°N 

A0 

V 

A 

0 

E 

'4> 

U,T,n) 


=  r  exp  i  ( 4>  —  tt/2) 

coordinates  in  transformed  space  (eq.  (3)) 

angle  of  attack 

density 

ratio  of  specific  heats 
meridional  coordinate  (fig.  4) 
normalized  radial  distance,  (r  -  rj3)/(rs  ”  rjj) 
damping  coefficient,  eq.  (18) 
eigenvalues  of  matrix  M 
eigenvalues  of  matrix  N 
turning  angle  defined  by  eq.  (34) 
gradient  operator 
difference  operator 
■  ♦  "  */2 
small  increment 
wing  leading  edge  sweep  angle 
coordinates  in  computational  space,  eq.  (7) 


Subscripts 

r,m  radial,  meridional  direction 

1,2  pre  and  post  shock  conditions 


Superscripts 

dummy  parameters  defined  following  eq.  (18) 
T  transpose  operator 
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Wing-body  angle  of  attack,  a  ,  deg 

Figure  1.-  Effect  of  angle  of  attack  and  free-stream 
Mach  number  on  interference  of  body  on  wing. 


Section  A  -  A 
z  *  10.04 


Figure  2.-  Flow  in  crossflow  plane  of  body  of  revolution 
with  supersonic  crossflow  Mach  number  as  obtained 
from  Euler  equations. 


Bow  shock 


_J _ 

28.  5 

Figure  3a.-  Mesh  for  delta  wing  alone 


Figure  3b.-  Detail  of  mesh  near  wing  leading  edge 


Transformed  Space 

(2  ,r,*) 


Figure  4.-  The  transformation  between  the  physical  space 
and  the  arbitrary  curvilinear  space. 
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Shock 


(b)  Computational  plane. 


Figure  5.-  Continued. 
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Figure  6.-  Steady  shock  for  three-dimensional 
supersonic  flow  in  generalized  curvilinear 
coordinates  (r,r,<£). 


Aspect  ratio:  1.0 
Thickness:  0.5 

Root  chord:  12  in 
Bevel  half  angle:  15° 

7a.-  Pressure  wing  geometry. 


Figure  8.-  Conical  streamlines  for  delta  wing  alone 
computed  by  Euler  code  with  Kutta  condition. 


PARTICLE  PATHS 


VELOCITY  VECTORS 


x/S 


11a.-  Pressure  on  body  at  a  station  just 
prior  to  start  of  wing. 
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Leeside  pressure 
Windward  side  pressure 


Pressure  on  wing-body  at  the  trailing  edge. 


interference  factor 
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